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We study Q-ball formation in the expanding universe on ID, 2D and 3D lattice simulations. 
We obtain detailed Q-ball charge distributions, and find that the distribution is peaked at Qp^k — 
1.9 X lO^^d^inl/m)^, which is greater than the existing result by about 60%. Based on the numerical 
simulations, we discuss how the Q-ball formation proceeds. Also we make a comment on possible 
deviation of the charge distributions from what was conjectured in the past. 
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I. INTRODUCTION 

If scalar fields are ubiquitous in nature, some of them may remain light in the early universe. Such scalar fields 
could be deviated from the low-energy minimum due to quantum fluctuations during inflation and/or modification 
of the scalar potential through gravitationally suppressed interactions with the inflaton. At a later time when the 
expansion rate becomes smaller than the mass, the scalar will start to roll down and oscillate about the potential 
minimum. The dynamics of such a scalar can exhibit rich phenomena in cosmology. 

We focus on a complex scalar field with a global U(l) symmetry at the potential minimum. The U(l) symmetry is 
not necessarily exact, and could be violated at high energy scales; our arguments apply if it is a good symmetry at 
low energy. If supersymmetry (SUSY) is realized in nature, there may be many such scalar fields. In Ref. [1 Coleman 
found that, if the scalar potential is shallower than a quadratic potential, there is a non-trivial field configuration that 
minimizes the energy for a fixed U(l) charge Q. Since the solution has a spherical symmetry, the solution is named 
Q-ball. A Q-ball is a non-topological soliton whose stability is supported by the U(l) symmetry. 

It was noticed in Refs. [2]-[5] that Q-balls play an important role in a context of the Affleck-Dine (AD) mechanism [6]. 
The mechanism utilizes a fiat direction of the supersymmetric standard model (SSM), which possesses a non-zero 
baryon (or lepton) number. A flat direction responsible for the AD mechanism is referred to as the AD field (denoted 
by $ in the following). The AD field ^ develops a large expectation value during inflation, and it starts to oscillate 
after inflation when the cosmic expansion rate becomes comparable to its mass. The baryon number is effectively 
created at the onset of the oscillations. Finally, <I> decays into the ordinary quarks, leaving the universe with a right 
amount of the baryon asymmetry. It was realized however that, soon after the onset of oscillations, the AD field 
experiences spatial instabilities and deforms into clumpy Q-balls [4] . Later it was shown in [7 that most of the baryon 
asymmetry is absorbed into Q-balls. Its non-linear property necessitates numerical approach to the formation and 
the subsequent evolution. 

The properties of Q-balls crucially depend on the size of the U(l) charge Q. For instance, the mass per unit charge is 
a decreasing function of Q in case of gauge-mediation type Q-balls; for a large enough Q, the Q-ball can be absolutely 
stable against the decay into nucleons and therefore contribute to dark matter [4J. In the case of unstable Q-balls, 
the decay rate depends on the charge, and Q-balls can be very long-lived in a cosmological time scale for a sufficiently 
large Q [8], since the Q-ball decay process takes place only in the vicinity of the surface [9 . Thus it is of utmost 
importance to determine the charge distribution of the Q-balls at the formation. 

There is a number of numerical studies on the Q-balls. The Q-ball formation has been first examined by Kasuya and 
one of the authors (MK) on 3D lattice simulations in a context of gauge mediation [TOl [11] and gravity mediation [7] , 
taking account of the cosmic expansion. They found that most of the baryon number is absorbed into the largest 
Q-balls in the lattices and estimated the typical Q-ball charge. However, the estimated charge might have been 
affected by relatively coarse spatial resolution due to the limited computational power at that time. The Q-ball 
charge distribution was studied in 2D ^ and 3D [13 lattices, with an emphasis on negative Q-balls formed for an 
orbit with a large ellipticity. In gravity mediation, however, the AD field naturally acquires a large kick into the 
phase space, leading to an orbit with a comparatively small ellipticity. As we will see later, the formation of negative 
Q-balls are also observed in our simulations for a large ellipticity. However, the final charge distribution does not fit 
well the empirical formula conjectured in Ref. [12 based on the entropy argument. The difference may be attributed 
to the fact that the charge re-distribution through exchanges of small secondary Q-balls seems less efficient in our 
simulations. Also the finite grid size could have affected the Q-ball evolution especially at a very late time when one 
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Q-ball was represented by a single grid point. In our analysis we stop following the evolution once the spatial resolution 
becomes insufficient to describe the Q-ball solution, in order to avoid any artificial effects due to the finite grid size. 
Recently Tsumagari studied the Q-ball formation and thermalization processes in a great detail in a Minkowski space 
neglecting the cosmic expansion [14 . In an expanding universe, however, the thermalization processes through Q-ball 
interactions may be weakened, and the resultant Q-ball configuration could be quite different. We will see in the 
following that the cosmic expansion indeed plays an essential role in the Q-ball formation. 

In this paper we study the Q-ball formation in a scalar potential motivated by the AD mechanism in gravity 
mediation and obtain the Q-ball charge distribution on ID, 2D and 3D lattice simulations with the cosmic expansion 
for different values of ellipticity and strength of instabilities. We develop a sophisticated numerical code based on the 
6th-order symplectic integrator equipped with a new Q-ball identification algorithm. 

The paper is organized as follows. In Sec. |TI| we summarize the governing equation of Q-balls in the gravity- 



mediation model. In Sec. Ill, we briefly mention the numerical scheme used in the simulations and, in Sec. |IV[ we 



show numerical results obtained from simulations in ID, 2D and 3D configurations. After that, Sec.[V|is devoted to 



some applications of our numerical results, and we conclude in Sec. VI 



II. Q-BALL DYNAMICS 



In the gravity-mediated SUSY breaking model, the AD field typically has a potential of the following form. 



1 + i^log 



1^ 



2n-2 



(1) 



where $ is a complex scalar field parametrizing the flat direction, A a coupling constant of the non-renormalizable 
term, H the Hubble parameter, c a positive constant, M a cutoff scale of the theory, the renormalization point 
defining the mass of scalar field m, and K = —0.01 ~ —0.1 is the contribution from the one- loop correction chiefly 
due to gauginos [5 . The cut-off scale M is typically taken as the GUT scale ~ 10^^ GeV or the Planck scale. Note 
that the potential ([T]) is invariant under the phase rotation of <I>. The effect of the U(l)-violating term will be taken 
into account in the initial condition of <I>. 

The evolution equation of the scalar field <I> is given by V^V^$ = —dV{^)/d^^. With the above potential, the 
explicit form of the governing equation of <I> becomes 



(n - 1)A^ 



^$ = 0, 



(2) 



where D is the spatial dimension, the dot indicates a time derivative, and is defined as the Laplacian with respect 
to the comoving coordinate. In this paper, we restrict our analysis to the epoch between the end of inflation and 
the reheating era, where the scale factor behaves as a cx t^/^ like the matter-dominant case. In addition, the non- 
renormalizable term (the last term being proportional to A) does not directly contribute to the late-time evolution of 
the scalar field, so we fix A = 0, and, for simplicity, c = 1. 
The charge density of the AD field is given by 



q{t,x) = -z(<l>*^-^$*). 
The charge and energy of a Q-ball are given by 

q{x) d^x^ 



E 



-L 



|$|2 + |V$|2 + y($) 



(3) 

(4) 
(5) 



where V denotes the /^-dimensional volume of the Q-ball. Note that, from the dimensional analysis, these quantities 
in ID and 2D should be regarded as the densities per unit area and length, respectively. 

For a negative value of the potential ([l]) is shallower than the quadratic potential m^|<l>p, for which there exists 
a non-trivial field configuration that minimizes the energy for a fixed charge Q. We can find the solution using the 
Lagrange multiplier method, and the minimum energy state turns out to be a spherical field configuration inside 
which the scalar field is rotating in a circular orbit in the phase space [1]. This object is called Q-ball. 
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There is an instability band for the potential with K < 0. From a simple linear analysis [7] one can show that the 
instability band is given by ^ 

< ^ < 2m^\K\. (6) 

The small fluctuations within the above instability band will therefore grow into Q-balls with a typical size 
R ~ |i^|~^/^m~^. Since the Q-ball formation is a highly nonlinear phenomenon, we need to rely on numerical 
simulations to understand how it proceeds. 



III. SIMULATION SETUP 

The potential ([T]) is invariant under the U(l) rotation of and therefore the charge Q is preserved. In general, 
however, there is a U(l)-violating term arising from the A-term oc {^'^ + <|)*^). This term kicks the AD field into the 
angular direction at the onset of oscillations, giving a nonzero U(l) charge, or equivalently, the angular momentum in 
the complex <l> plane [see Eq. (|3|]. Since the amplitude of the oscillations decreases due to the Hubble friction term, 
such non-renormalizable U(l)-violating terms soon become negligible. Thus, the main role of the A-term is to give 
an initial angular momentum. To take account of the effect, we parametrize the initial condition of the homogeneous 
mode as 

$(tin) = M*, ^(tin) = ieM^m. (7) 

The initial charge density is then q{tin) — 2eM^m. In the following analysis we take e = 0.01, 0.1, and 1. Also we add 
small fluctuations to the above initial condition, |(5<l>/<l>| = 10~^, following Ref. [7 . 

We use a comoving box as a computational domain whose physical size at the initial time is a(t\n)L = hH^ ^ where 
i^in is the Hubble parameter at the initial time, and 6 is a dimension-less parameter defining the box size. We choose 
the normalization of the scale factor as a{t\n) = 1 where is the initial time, = 2/{?>H\jX We adopt a value 
between 6 = 1 ~ 5. In all simulations except for the partial cases (1D2, 1D3, 1D4 in Table m in ID, we use the 
equation of motion in an expanding background with an initial condition H-^^ = m. In the non-expansion case, L 
coincides with the physical scale. As for the computational grid, we use the homogeneous N^^-^^ grid, where D is the 
dimension and 256 < A/gHd ^ 4096 for D = 1, A/grid = 512 for I) = 2 and A^gnd = 128 for = 3. Other numerical 
parameters we used are tabulated in Tables IpII On the boundaries, we impose the periodic boundary condition. 



As the universe expands, the spatial resolution gets worse. The physical Q-ball size is given 
and the grid spacing is Ax = a{t)L/Ngrid' Requiring R > Ax during the simulations, we obtain the critical time after 
which Q-balls are no longer resolved, 

.,.mfe-U = ^(j^y'\ (8) 

where r measures the time normalized by from the beginning of the simulation, i.e., r = m(t — tin). For example, 
we can resolve Q-balls until r < 1.1 x 10^ for A^grid = 1024 and 6 = 5. In all simulations we have chosen the number 
of grids so that Q-balls are resolved until r ~ 10^. 

To evolve <I>, we implement the 7 stages 6th-order symplectic integrator developed by Yoshida [16]. This scheme 
conserves the charge and energy of the target field, and the higher-order scheme is expected to suppress the global 
error, at least, until the epoch of the Q-ball formation. The spatial derivatives, i.e. the Laplacian operator in Eq. ([2|, 
are estimated with the second-order finite difference. We also implemented the Fourier collocation method, and 
confirmed that the final results do not depend on the scheme of derivative. 

In order to identify the formed Q-balls, we developed an algorithm based on the Marching Cube method [17]. 
We briefly summarized the algorithm in Appendix B. Our algorithm regards the region in which the charge density 
exceeds a critical value, \q{t^ x)\ > as a Q-ball. This criterion is determined so that Qc is sufficiently larger than the 
background field but smaller than the typical peak density. We set Qc to be a fraction of q{t) when the fluctuations 
first enter the non-linear regime. 

In addition, we use OpenDX to visualize the Q-balls in the 2D and 3D cases [18]. 



In Ref. [7], there is an error in the instabihty band. One of the authors has rederived the instabihty band with a more plausible assumption 
in Ref. [15]. 
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IV. NUMERICAL RESULTS 
A. ID 

1. the effect of cosmic expansion 

Let us begin with the time-evolution of the AD field in ID. For the sake of seeing the effect of the cosmic expansion 
particularly in ID, we have performed simulations both with and without the cosmic expansion. 

First, we show the result with the cosmic expansion and explain how to identify Q-balls. Fig. [l] shows the charge 
density of the AD field at r = 800 and 5000 with the parameter set IDl (see Table There are lots of peaks to be 
regarded as Q-balls. These Q-balls were formed at r ~ 700 in this simulation when the background charge density 
has been decreased to gform ~ go/<^(700) ~ 0.02 \ ^in\^m where is the initial value of ^ given in Eq. Once a 
Q-ball is formed, its charge density at the center remains virtually unchanged unless there are significant collisions 
and mergers, since its physical size does not change in time and the total charge is preserved under the theory we 
consider here. Actually, we can see that the peak values are more or less unchanged between r = 800 (red solid) and 
r = 5000 (green dashed) in Fig. [l] 

At the formation, the modes in the instability band ([6| grow and become non-linear, and proto-Q-balls are formed. 
Those proto-Q-balls absorb the surrounding charges, and the charge density outside Q-balls are quickly damped. 
Therefore if we set qc to be a fraction of g'form, we can collect the Q-balls. Throughout this paper, we choose 

= ^form/5. The explicit value of qc depends on each simulation. In the present case, = 4 x 10~^|<l>inpm, which 
is drawn as the blue dotted line in Fig. [T] 

Using this criterion, we show the time evolution of the number of Q-balls in Fig. [2] The vertical axis is the number 
of Q-balls in unit (comoving) length of the computational domain, A^tot- We found that there are approximately 5 Q- 
balls in the unit comoving length (normalized by m), that is, totally ~ 25 Q-balls in the whole domain with 6 = 5. To 
reduce the statistical error, we averaged 50 realizations. We can see that, after Q-balls are formed at about r ~ 700, 
the number of Q-balls remains constant for a while and then gradually decreases. The decrease is probably due to 
the merger or disruption process through collisions. In the ID simulations, collisions among Q-balls are expected to 
occur more frequently than in the cases with 2D or 3D, which enhances the probability that multiple Q-balls merge 
into a large one. Note that this is not due to the finite box effect, as will become clear from the following argument. 



ID 


# of grid 


Lm 


dt 


K 


expansion 


realization 


IDl 


1024 


5 


0.005 


-0.1 


yes 


50 


1D2 


256 


100 


0.05 


-0.1 


no 


50 


1D3 


1024 


400 


0.05 


-0.1 


no 


50 


1D4 


4096 


1600 


0.05 


-0.1 


no 


20 



TABLE L Numerical parameters for ID simulations. 



Identifying each Q-ball, we can estimate its charge and energy. The charge distribution of Q-balls at r = 1000 
(blue dotted), r = 5000 (green dashed) and r = 8000 (red solid) is shown in Fig. [s] The vertical axis is the number 
multiplied by the charge. Note that the number here means the one in a logarithmic bin in a unit comoving length, 
N{Q). The charge distribution does not significantly change in time except for a temporal bimodal feature at r = 1000, 
and the Q-balls mainly contributing to the total charge in the computational domain are those with Q :^ l^inP- Note 
that Q should be interpreted as the charge density per unit area, and so, it has a mass dimension two. Moreover we 
have found that the Q-ball charge ranges over about one order of magnitude, 0.2 < Q/|^inP ^ 2. 

In Ref. [7 , it was concluded that based on the 3D simulation, a small number of very large Q-balls absorb most of 
the charge that the AD field originally possessed. We however have found that the charge has a distribution with a 
finite width shown in Fig. [3j Therefore, in this paper, we regard the typical charge of Q-balls as that at the peak of 
the charge distribution NQ instead of the rare largest Q-balls. 

So far, we consider the case with the cosmic expansion. In the non-expansion case, the spatial resolution does not 
change during simulations. So, it may sound good for the purpose of numerical simulation and there are some past 
works assuming the non-expanding background, e.g. [14]. However, there are several crucial differences. 

Fig. [4] shows the time evolution of the number of Q-balls in the non-expansion cases with the parameters 1D2, 1D3 
and 1D4 tabulated in Table |lj The red solid, green dashed and blue dotted lines are the results with a different box 
size: Lm = 100, 400 and 1600, respectively. The vertical axis is the total number of Q-balls in a simulation box. 
These simulations have the same spatial resolution Ax = Lm/N^^i(\. In order to reduce the statistical fiuctuations, we 
average over 50 realizations for Lm = 100, 400 and 20 realizations for Lm = 1600. The Q-balls are formed at r ~ 400, 
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0.15 




0.05 



FIG. 1: The spatial charge-density distribution in the simulation IDl. The blue dotted line is the criterion we set for this 
simulation, given by = 4 x 10~^\^in\^m. The horizontal axis is the comoving coordinate. 
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FIG. 2: The time evolution of the total number of Q-balls in the unit comoving length, for the parameter set IDl and 



and their number monotonically decreases due to merger and disruption through collisions. The decrease may be 
similar to what we observed in the case with the cosmic expansion. However, at about r-2xl0'^,7x 10^,1x10^ for 
Lm = 100, 400, 1600, respectively, the Q-ball number starts to rapidly decay. In fact, the numbers tend to converge 
to 1, and the convergence time as well as the time when the rapid decay starts get delayed for a larger box. This 
observation shows that the large Q-balls sweep the whole computational domain many times, absorbing other small 
Q-balls, until one very large Q-ball is formed. This is clearly the artificial volume effect. Therefore the numerical 
results in the non-expansion cases can be trusted only before the time when the numbers start to rapidly decay. 

The difference between the non-expansion and expansion cases can be manifestly seen from Fig. |5j This figure 
shows the trajectories of Q-balls in the case with the cosmic expansion (IDl in Table |l| [right] and without (1D3) 
[left]. In the non-expansion case, we found that the Q-balls interact many times with one another, and Q-balls go 
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FIG. 3: The charge distributions of Q-balls at r = 1000 (blue dotted), r = 5000 (green dashed) and r = 8000 (red solid), for 



the parameter set IDl and = 4 x 10 
length, N{Q), multiplied by the charge. 



|$in| m. The vertical axis is the number of Q-balls with charge Q in a comoving 




FIG. 4: The evolution of the total number of Q-balls in the unit comoving length in non-expansion cases. The box size Lm 
are 100, 400 and 1600 while the spatial resolution of these simulations are fixed. For r > 10^, we expand the output interval to 
avoid the dense plot. 



across the boundaries of the domain again and again. Through the interactions among Q-balls, large Q-balls tend to 
become even larger by ingesting other slow-moving Q-balls, resulting in the decrease of the number of Q-balls. Since 
the background spacetime is static, the process will never stop until all the charge is absorbed into one large Q-ball. 
Hence, in order to avoid the artificial volume effect, a large computational box must be prepared. 

In the expanding background, Q-balls are dragged by the cosmic expansion and thus less frequently interact with 
other Q-balls. In fact, as shown in the right panel of Fig. |5j the number does not change so much even with a quite 
small box. The cosmic expansion clearly suppresses the interactions among Q-balls, and the merger and disruption 
processes do not frequently take place; they are considered to be decoupled at a certain point. In the 2D and 3D 
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FIG. 5: Trajectories of the Q-balls in the cases without (left) and with (right) the cosmic expansion. Each dot represents the 
position of a Q-bah. The horizontal axis of the right panel is the comoving coordinate, and that of the left one is the physical 
coordinate in unit of . 



cases, the finite box effect may look weaker than in the case of ID, since the collisions are expected to take place less 
frequently. However, the intrinsic problem of neglecting the cosmic expansion mentioned above still exists even in 
these cases, and it is quite likely that such a crude approximation significantly affects the Q-ball formation and the 
final Q-ball charge distribution. Thus we hereafter will treat only the cases with the cosmic expansion. 



B. 2D 



1. Filamentary structure 

Let us next consider the 2D system. The parameter sets used in our analysis are shown in Table [Til We first show 
the result with the parameter set 2D1. The snapshots of the spatial charge distribution during the Q-ball formation 
are shown in Fig. 6j In this figure, the panels (a)-(d) correspond to r = 725,825,875 and 2000, respectively, and the 
red indicates the large positive charge and the black q ^ 0. We can clearly see the filamentary structure just before 
the Q-ball formation [see (a)(b)], which was noted in Ref. [12]. Our point here is that the filamentary structure plays 
an important role in the Q-ball formation. In order to help understanding the formation process, we classify the 
Q-balls into those formed at intersections, sides and voids of the filamentary structure, in the order of the charges 
and the formation time. ^ At r ~ 800, large proto- Q-balls are formed at the intersections of the filaments, as one 
can see in (b). These proto-Q-balls will chiefly become the largest Q-balls forming the peak in the final distribution. 
Shortly after the formation of the proto-Q-balls at the intersections, the filaments start to be torn to smaller pieces, 
which will in the end become Q-balls with smaller charges. In addition, small Q-balls are formed in the void regions 
of the filamentary structure almost at the same time [see (c)]. The Q-balls formed from the sides and the voids will 
account for the tail of the charge distribution left to the peak. These classification may become vague after mergers 
and disruptions through collisions, but it is clear from Fig.|6]that the filamentary structure strongly affects the spatial 
charge distribution even at late times. This argument is reinforced by the observation that major collisions are less 



^ Note that this classification is obscure especially for small Q-balls; significant Q-ball interactions are considered to occur locally, which 
makes it difficult to chase the origin of those Q-balls. However it is still useful to keep in mind such classification because those small 
Q-balls do not significantly contribute to the peak of the charge distribution. 
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frequent in the expanding background, as we have seen in the case of ID (see Fig. |5|. Also we can exphcitly see that 
the number of Q-bahs is preserved soon after the Q-bahs from the sides and voids are formed, as shown in Fig. [7| 



ID 


grid size 


b dt 


K 


e 


2D1 


512^ 


2.5 0.005 


-0.1 


1 


2D2 


512^ 


2.5 0.005 


-0.07 


1 


2D3 


512^ 


3.0 0.005 


-0.04 


1 


2D4 


512^ 


2.5 0.005 


-0.1 


0.1 


2D5 


512^ 


2.5 0.005 


-0.1 


0.01 



TABLE II: Numerical parameters for 2D simulations. 




(c) r = 875 (d) r = 2000 

FIG. 6: Time evolution of charge density of the AD field in 2D. The red indicates the large charge region. For these visualizations, 
we used OpenDX. 



2. Charge distribution 

Fig. [s] shows the charge distribution of Q-balls for the 2D1 simulation. We set the criterion = 3 x 10~^|<l>inpm 
for the identification of Q-balls. The distribution does not significantly change once the Q-balls of the three types 
are formed. This is consistent with our observation that major collisions are rare in the expanding universe. The 
distribution has a sharp peak at Qp^^k ^ 0.1|^inP^~^ and the typical charge ranges over one order of magnitude 
0.03 <Q/(|$inpm-i)< 0.2. 
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1000 



FIG. 7: The time evolution of the total number of Q-balls in the unit comoving area, for the parameter set 2D1 and Qc 
3 X 10-^|$in|^m. 



Let us compare the existing result with ours. The maximal charge observed in the simulation of Ref. [TT] reads 

Qmax"^ ~ 0.1|<l>inp7Ti~^ in our sctup^. This agrees very well with the above peak value. 

We confirmed that the Q-balls formed at the intersections account for the largest Q-balls, especially those around 
the peak, while the Q-balls created from the sides and the voids form a long tail toward smaller charges. Thus we 
conclude that it is the filamentary structure that determines the charge distribution, for an orbit with a small ellipticity. 



0.6 




FIG. 8: Charge distributions of Q-balls at r = 2000 (green dashed) and r = 5000 (red solid), for the parameter set 2D1 and 
= 3 X 10~^|$in|^m. The vertical axis is the number of Q-balls with charge Q in a comoving area, N(Q), multiplied by the 
charge. 



^ Please notice that in Ref. [J, the treatment of ID and 2D Q-balls is different from ours. The 2D Q-balls in Ref. [? ] has an infinitely- 
long sylindrical shape in 3D space. 
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3. Relationships among charge, energy and size 



The gravity- mediation type Q-ball satisfies a relation between the charge and the energy, E ~ mQ. To check 
this relation, we plot in the left panel of Fig. [9] the charge and the energy of the Q-balls in the 2D1 simulation. 
This figure shows that these quantities are strongly correlated, having an almost linear relation, |(5|/(|<^>inP^~^) — 
0.42(E/|$inp)°-^^ for IQI > 10-^|cE>inpm-\ as expected. Thus, it confirms that our numerical algorithm for the 
formation and the identification of Q-balls work well. We note however that the proportional coefficient is smaller 
than 1, which means that the orbit of the AD field in the phase space is not completely circular. Strictly speaking, 
those objects we have identified as Q-balls are not the Q-ball solution but slightly excited states of the Q-balls [19]. 
This excitation may be due to the Hubble friction at the early time for the simulation. At the initial time, we set the 
initial field so that mQ/E ^ 1 [see Eq. The field, however, cannot generally keep the circularity of its orbit in the 
phase space while the Hubble friction works efficiently. Also, we will see the dimensional effect on the proportional 
coefficient in Sec. IIV C 31 

We have seen no further relaxation process toward the true Q-ball solution in our numerical simulations. This 
is a limitation of our analysis. However, since the coefficient is relatively close to 1, we expect that those objects 
will gradually approach the Q-ball solution in a much longer time scale ^, without significantly affecting the charge 
distribution obtained here. 

Before showing the relation between the charges and the sizes, we mention the definition of diameter of a Q-ball in 
this paper. The shape of a Q-ball identified by our algorithm is generically a distorted disk (ball in the 3D system). 
So, our procedure to calculate the diameter oi a Q-ball is as follows. First, we measure the area (the volume in the 3D 
system) of the region on the grid in which the charge density is larger than the criterion. At this stage, the location 
of the boundary of a Q-ball is specified by the linear interpolation using the two values of q{t^ x) at the adjacent grid 
points between which q(t^ x) = is realized. Then the definition of the diameter of a Q-ball in this paper is the 
diameter of the disk (a ball) with the same area (volume) as measured above, which is given in Eq. (Bl). For details, 
see Appendix B. 
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* positive 
X negative 

fitting : 0.42E' 

x=5000 




FIG. 9: The charge and the size (diameter d = 2R) of Q-balls with parameter set 2D1 and = 3 x 10~^\^in\^m. The red 
'+' and the blue 'x' represent positive Q-balls and negative ones at r = 5000, respectively. The green dashed line represents 
the fitting formula, |Q|/(|$in|^m~^) = 0.42(£^/|$in|^)°"^^. The arrow standing on the horizontal axis represents the grid size 
at r = 5000 (Ax 1.87m"^). 

In the right panel of Fig. [9j we plot the relation between the size (diameter d = 2R) and the charge. The Q-ball 
solution has a certain size determined from the instability band shown in Eq. (|6|. We have found that all the Q-balls 
at least above Q > 10~^|<l>inpm~^ have an almost same size, d ~ 10m~^. This is another evidence that our numerical 



The relaxation time may be extremely long due to an approximate adiabatic invariant |2Q] . 
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scheme correctly follows the Q-ball formation and collects the Q-balls. The spatial resolution in this simulation is 
presented as an arrow standing on the horizontal axis of Fig. |9j being much smaller than the typical size of a Q-ball. 
Hence the Q-balls are successfully identified on the lattices with sufficient spatial resolution. 



4- K- dependence 

The strength of the instability is controlled by the coefficient of the one-loop corrections, K. For a larger value of 
\K\^ the instability gets stronger and the Q-ball size R ~ |i^|~^/^m~^ becomes smaller. Thus, the final Q-ball charge 
distribution may depend on the value of K. 



Figure [10] shows the charge distributions at r = 5000 with K = —0.1, —0.07 and —0.04, corresponding to 2D1, 2D2 
and 2D3, respectively, in Table [ll} Although the charge distribution seem to shift toward larger values of Q as 
no significant difference was observed. 
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K=-0.04 
K=-0.07 
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T=5000 
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FIG. 10: Charge distributions of Q-balls with i^^ = -0.1, = 3x 10-^|$in|^m [2D1] (red solid), K = -0.07,^c = 2xl0-^|$i, 
[2D2] (green dashed), and K = -0.04, = 1.2 x 10"^|$in|^m [2D3] (blue dotted). 



The size (diameter d = 2R) and the charge of Q-balls in simulations 2D1-2D3 are shown in Fig. 11 Focusing on the 
large Q-balls with Q = Qp^^k ^ -^^^ each result ^, we can read the typical sizes as (i = 2R ~ 14m~^, 17m~^, 24m~^ 
for K = —0.1, —0.07, —0.04, respectively. Thus we confirmed that the size approximately scales as ~ |i^|~^/^m~^. 
More specifically, we obtain R 2.3|K|"^/^m"^ at r = 5000. 



5. Initial elliptic orbit 

So far we have concentrated on the cases with e = 1 [see Eq. ([t]) for its definition] . If the parameter e is smaller than 
1, the Q-ball formation and the charge distribution are significantly changed. However, the filamentary structure is 
still essential to understand the final distribution. 

For e ^ 1, the initial orbit of the AD field becomes elliptic. It has been known that negative Q-balls, namely, 
Q-balls with a negative charge, appear in this case [7l|12]. This is because the AD field cannot grow into the Q-balls 
while keeping the energy-to-charge ratio E/mQ ~ 1/e ^ 1. (Recall that the Q-ball solution satisfies E ~ mQ.) 
Thus, the AD field must somehow discard its excessive energy to reach the Q-ball solution, while the total charge 



^ The reason why we select the large ones is that the error in measuring the size is considered to be minimized. Smaller Q-balls tend to 
have smaller radius in our Q-ball identification algorithm due to the presence of the ambient background fiuctuations. 
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FIG. 11: The charge and the size (diameter d = 2R) of Q-balls with K = -0.1, = 3 x 10"^|$in|^m [2D1] (red '+'), 
K = -0.07,^c = 2 X 10"^|$in|^m [2D2] (green 'x'), and K = -0.04,^c = 1.2 x 10"^|$in|^m [2D3] (blue The arrows 
standing on the horizontal axis represent the grid size at r = 5000 for K = —0.1,0.07 simulations (tall one. Ax = 1.87m~^) 
and for K = —0.04 simulation (short one. Ax = 2.25m~^). 



is conserved. This can only be achieved by changing the ratio locally; in reality, this results in producing negative 
Q-balls. 

We have performed numerical simulations with e = 0.1 (2D4) and e = 0.01 (2D5). The snapshots of the spatial 
charge distribution at r = 800, 1500,2200 and 5000 are shown in Figs. [12] and [Tsj respectively. In these figures, the 
red represents the positive Q-balls and the cyan the negative ones. Let us first explain how the formation proceeds 
in the case of e = 0.1. As in the case of e = 1, the filamentary structure appears first [see the panel (a)]. As the 
field grows, the filaments are torn to smaller pieces and the proto-Q-balls are formed from intersections, sides and 
voids of the filaments. Similarly to the case of e = 1, the largest Q-balls are those from the intersections. At this 
stage, we call these Q-balls first- generation Q-balls. Almost all the first-generation Q-balls have a positive charge, 
as can be seen in the panel (b). Note that the first-generation Q-balls are in a highly excited state in a sense that 
they have large E/Q ratio. Around r ^ 2200 [the panel (c)], negative Q-balls (cyan) begin to be formed. The point 
to be emphasized here is that the negative (as well as positive) Q-balls, especially larger ones, appear around the 
existing first-generation Q-balls. If one takes a closer look, one may notice that a pair of positive and negative Q-balls 
is produced from one first-generation Q-ball. This observation clearly shows that the highly excited first-generation 
Q-balls release the excessive energy in a form of positive and negative Q-ball pairs. At r = 5000 [the panel (d)], the 
relaxation process has been almost finished. We call these Q-balls second- generation Q-balls to distinguish them from 
the first-generation Q-balls. We can see some isolated negative second-generation Q-balls, but most of them are still 
paired with the positive ones. 

In the case of e = 0.01, the evolution is quite similar until the first-generation Q-balls are formed (the panel (b) 
in Fig. 13). The production of the negative Q-balls are more violent as can be seen from the panel (c). We can see 
multiple positive and negative Q-balls are produced from one first-generation Q-ball. This should be contrasted to the 
previous case that a pair of positive and negative Q-balls is produced. Most of the large negative and positive Q-balls 
remain paired even at r = 5000. Note that the largest (positive or negative) Q-balls are produced from the largest 
first-generation Q-balls, which chiefly arise at the intersections of the filamentary structure. In the void regions, many 
small Q-balls appear in comparison with the case of e = 0.1. This may be induced by the energy released from the 
first-generation Q-balls in the violent relaxation process. 

We expect that the filamentary structure also affects the final charge distribution, since it played an essential role 
in the Q-ball formation. Fig. [14] shows the charge distributions of positive and negative Q-balls. The upper panels are 
the result with e = 0.1, and the lower ones that with e = 0.01. We take the critical charge Qc being proportional to e 
as Qc = ^3 X 10~^e|^inpm, since the initial charge density is proportional to e. In the case of e = 0.1, most Q-balls 
are positive, and negative Q-balls are subdominant. The peak charge is approximately Qp^ak ^ 0.015|$inpm~^ which 
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(c) T = 2200 (d) T = 5000 



FIG. 12: The color map of the charge density for e = 0.1. The red represents the positive charge, the cyan the negative one 
and the black ^ ^ 0. 

indicates the scaling as Qp^^k ^ ^- (R^ecall that the peak charge is Qp^^k ^ 0.1|^inP^~^ for e = 1. See Fig. [sj) On 
the other hand, in the case of e = 0.01 where the initial orbit is highly eccentric, almost the same number of positive 
and negative Q-balls are produced after r ^ 2000. The peak charge is approximately 2 x 10~^ until r < 2000 which 
agrees with the scaling. However, the scaling of the peak charge is no longer valid at r = 5000 [TT]. 

In order to study how the Q-balls in an excited state approach toward the Q-ball solution, we show the energy 
and the charge of the Q-balls in the left panels of Fig. [Ts] and [l6j and the relation between the size (diameter) 
and the charge in the right panels. In the figures, the green symbols ('*') represent the energy and charge of the 
first-generation Q-balls at r = 1500, implying E/mQ 1/e. Besides, the dependence of the size on the charge is 
somewhat large. These facts indicate the first-generation Q-balls shown here are highly excited states from the static 
Q-ball solution. The red symbols ('+') and the blue ones ('x') represent the positive and the negative Q-balls at 
r = 5000, respectively. We can see that the energy-to-charge ratio decreases as time goes, while roughly keeping the 
proportionality between E and Q, and the sizes become almost universal. The Q-balls seem to stop evolving after 
the ratio E/mQ reaches 0(10). These figures manifestly show the relaxation process from the Q-balls in an excited 
state to a state with a lower value of E/Q. 

The relaxation process of an excited state of Q-balls is a highly non-linear phenomenon, which can be studied only 
by numerical simulations. From our results shown above, we have reached conclusion that Q-balls in an excited state, 
namely having a large E/mQ ^ 1, continue to transform the excessive energy into multiple (positive and negative) 
Q-balls, until the E/mQ ratio becomes as small as 0(10)^. In the case of E/mQ ^ 10, or equivalently, e <C 0.1, 



The Q-balls might have some amount of the angular momentum. However, we looked into the kinetic, potential and gradient energy of 
the Q-balls and found that the angular momentum is not the main source for the large E/mQ ratio. 
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(c) T = 2200 (d) T = 5000 



FIG. 13: The color map of the charge density for e = 0.01. The red represents the positive charge, the cyan the negative one 
and the black ^ ^ 0. 

therefore, both positive and negative Q-balls are produced in the relaxation process, forming an almost symmetric 
charge distribution. 

In the case with e = 1, we pointed out that the ratio E/mQ is larger than 1 due to the Hubble friction. In the 
present case, however, the Hubble friction is no longer relevant at the formation of the second-generation Q-balls. 
Hence the fact that the second-generation Q-balls have somewhat large E/mQ may be due to the highly non-linear 
interaction during the secondary formation process. To find out the reason, more detailed investigations on the 
secondary formation process are required, which is beyond the scope of this paper. 

C. 3D 



ID 


# of grid 


h dt 


K 


e 


3D1 


128^ 


1.0 0.005 


-0.1 


1.0 


3D2 


128^ 


1.0 0.005 


-0.07 


1.0 


3D3 


128^ 


1.6 0.005 


-0.04 


1.0 


3D4 


128^ 


1.0 0.005 


-0.1 


0.1 


3D5 


128^ 


1.0 0.005 


-0.1 


0.01 



TABLE III: Numerical parameters for 3D simulations. 
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FIG. 14: Charge distributions of Q-balls with qc — ±3 x 10~^e|$in|^m. The upper panel is the case with e = 0.1 [2D4], and 
the lower one is e = 0.01 [2D5]. The left panels are distributions of negative Q-balls and the right panels are those of positive 
Q-balls. 



1. Filamentary structure 



In this final section, we move on to 3D simulations. The parameter sets used are tabulated in Table III 
The snapshots of the spatial charge distribution with the parameter set 3D1 are shown in Fig. [T7| The surfaces in 
the figures represent the isodensity surfaces of the charge density q = 2 x 10~^|<l>inpm (white), q = 1 x 10~^|<l>inpm 
(blue), and q = 5x 10"^|^inpm (red). The panels (a)-(d) correspond to r = 925,950,975 and 2000, respectively. As 
we have seen in the 2D results, we observed the filamentary structure that appears just before the Q-ball formation, 
and at each intersection, proto-Q-balls are formed represented by blue surfaces in the panel (a). At r = 950, Q-balls 
appear also on the sides, and at the same time, the Q-balls at the intersections grow which are represented by red 
surfaces. In the next panel (c), the filaments start to be torn to small pieces which eventually become small Q-balls. 
This panel shows manifestly that the Q-balls formed at the intersections hold large charges also in the 3D system, 
and thus they mainly contribute to the charge distribution. At the end of the simulation, r = 5000, we observed that 
the number of Q-balls in a unit volume becomes N = 10^ shown in Fig. 18 



The authors of Ref. [7] concluded that only a few number of Q-balls hold almost all charge of the whole simulation 
box. In contrast to this, our result indicates that, although the box size is smaller than that used in their simulations, 
there are quite large number of Q-balls with almost same size. This discrepancy may come from the difference of the 
spatial resolutions. In fact, while the resolution of the smaller box in their simulation is given by Ax ~ 0.05m~^ at 
the initial time, the one in our simulation is Ax ^ 0.008m~^. Moreover, even at the formation epoch, the resolution 
is given by Ax ~ a(lOOO) x m~^/128 ~ 1.0m~^, and thus Q-balls are resolvable at the epoch. In compensation for 
the fine resolution, the box size in our simulation is smaller than the past one. Nevertheless, Fig. [T7| indicates that 
the Q-ball formation caused by the collapse of the filamentary structure proceeds without being affected by the finite 
volume effect. 
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FIG. 15: The charge and the size (diameter d = 2R) of Q-balls with e = 0.1, gc = 3 x 10"^|$in|^m [2D4]. The green 
represents those of positive Q-balls at r = 1500. The red '+' and the blue 'x' represent positive Q-balls and negative ones at 
T = 5000, respectively. The arrows standing on the horizontal axis represent the grid size at r = 5000 (Ax = 1.87m~^). 
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FIG. 16: The charge and the size (diameter d = 2R) of Q-balls with e = 0.01, = 3 x 10~'^|$in|^m [2D5]. The green '*' 
represents those of positive Q-balls at r = 1500. The red '+' and the blue 'x' represent positive Q-balls and negative ones at 
T = 5000, respectively. The arrows standing on the horizontal axis represent the grid size at r = 5000 (Ax = 1.87m~^). 



2. Charge distribution 

Fig. [Tol shows the charge distribution of Q-balls for the 3D1 simulation. We set the criterion, Qc = 1.8 x 10~^|<l>inpm. 
Basically, the Q-ball charges have a distribution similar to the 2D cases, and those mainly contributing to the total 



charge range over an order of magnitude, namely, 2 x 10~^ < Q/(|^inP^~^) ^ 5 x 10~^. In Sec. VA, we will perform 
the fitting of this distribution in consideration of the Poisson error of the number of Q-balls in each bin. As a result, 
the typical Q-ball charge is Qp^^k — ^-^ ^ 10~^|*^inP^~^5 which is somewhat larger than the previous result obtained 
by Kasuya et al. [7 where they have concluded Qmax"^ = 1.2 x 10~^|<l>inpm~^ for K = —0.01, when expressed in 
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FIG. 17: The isosurfaces of charge density with = 2 x 10 ^\^u 



5 X 10"1$h 



(red). For these visuaUzations, we used OpenDX. 



(white), Qc 



1 X 10-^|$h 



(blue), and qc 



terms of our terminology. 



3. Relationships among charge, energy and size 



Following the previous discussion in Sec. |IVB 3l we investigate the relations among the charge and energy, and size. 
The left panel of Fig. [2Q| shows the charge and the energy of the Q-balls in the 3D1 simulation. Also in this 3D system, 
we can confirm that these quantities are strongly correlated, having an almost linear relation, \Q\/{\^i^\'^m~'^) 
0.17[E/{\^in\'^m~'^)]^-^'^ for \Q\ > 10~^|<l>inpm~^, as expected. Note however that the proportional coefficient is 
smaller than the one in the 2D cases. This small coefficient arises from the fact that the cosmic friction term in 
Eq. (§ is given by the spatial dimension. Thus the AD field tends to have a larger ellipticity in the 3D case compared 
to the 2D case, for a given initial condition in the phase space. As pointed out in Sec. |IV B 3[ these objects we 
identified as Q-balls may be slightly excited states of the Q-balls, and it is expected that those objects will gradually 
approach the Q-ball solution in a much longer time scale than the current limitation of our numerical analysis. 

In the right panel of Fig. 20, we plot the relation between the size (diameter d = 2R) and the charge (for the 
definition of the diameter, see the second paragraph in Sec. IV B 3). The Q-ball solution has a certain size determined 
from the instability band shown in Eq. (|6|. We have found that all the Q-balls at least above Q > 10~^|^inpm~^ 
have an almost same size, d = 2R ~ 10m~^, and that the size is insensitive to the dimension since this is almost 
same as the one in the 2D system. The spatial resolution in this simulation is presented as an arrow standing on the 
horizontal axis of Fig. [20| being much smaller than the typical size of a Q-ball. Hence the Q-balls can be successfully 
identified on lattices also in the 3D system. 
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FIG. 18: The time evolution of the total number of Q-balls in the unit comoving volume, for the parameter set 3D1 and 
qc = 1.8 X 10-^|$in|'m. 
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FIG. 19: Charge distributions of Q-balls at r = 2000 (green dashed) and r = 5000 (red solid), for the parameter set 3D1 and 
|$in|^m. The vertical axis is the number of Q-balls with charge Q in a comoving volume, N(Q), multiplied by 



Qc = 1.8 X 10" 
the charge. 



4. K- dependence 



The authors of Ref. [TT] have claimed that the largest charge of Q-balls in the 3D system depends on K as 



.(KK) 
^max 



(X 



\K\ 



-1/2 



To investigate this, we plot the charge distributions with K = 

Note that 



These parameter sets are shown as 3D1, 3D2 and 3D3 in Table III 



0.1,-0.07 and 
in case with K = - 
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-0.04 in Fig. 

-0.04, the size~of 

Q-ball becomes also larger. Hence we take a larger box than the other cases. 

From Fig. [2l] we cannot see any significant difference among different values of K. Suffice it to say that the 
distribution for K = —0.04 might be marginally skewed to large charges. However, to make this point clear, more 
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FIG. 20: The charge and the size (diameter d = 2R) of Q-balls with parameter set 3D1 and Qc = 1.8 x 10~^\^in\^m. The red 
'+' and the blue 'x' represent positive Q-balls and negative ones at r = 5000, respectively. The green dashed line represents 
the fitting formula, |Q|/(|$in|^m~^) = 0.17[£^/(|$in|^m~"^)]°"^^. The arrow standing on the horizontal axis represents the grid 
size at r = 5000 (Ax = 2.99m"^). 



realizations are required. If this is true, the deviation from the result of Ref. [T may rather widen because their result, 
QmJ^^ = 0.012|<|)inpm-2 for K = -0.01, would be smaller in the case with K = -0.1. 
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FIG. 21: Charge distributions of Q-balls with K = -0.1 [3Dl],^c = 1.8 x 10"'^|$in|^m (red solid), K = -0.07 [3D2],^c = 
1.2 X 10~'^|$in|^m (green dashed), and K = -0.04 [3D3],^c = 6 x 10~^|$in|^m (blue dotted). 



The size (diameter d = 2R) and the charge of Q-balls in simulations 3D1-3D3 are shown in Fig. 
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Focusing on the 
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large Q-balls with Q = Qp^^k ^ ^^^^ result, we can read the typical sizes diS d = 2R ^ 15m ^, 19m ^, 25m 

for = —0.1, —0.07, —0.04, respectively. Thus we confirmed that the size approximately scales as ~ |i^|~^/^m~ 
More specifically, we obtain R ~ 2.4|K|"^/^m"^ at r = 5000. 




FIG. 22: The charge and the size (diameter d = 2R) of Q-balls with K = -0.1 [3Dl],^c = 1.8 x 10"'^|$in|^m (red '+'), 
K = -0.07 [3D2],^c = 1.2 x 10-^|$in|^m (green 'x'), and K = -0.04 [3D3],^c = 6 x 10-^|$in|^m (blue '*'). The arrows 
standing on the horizontal axis represent the grid size at r = 5000 for K = —0.1,0.07 simulations (tall one, Ax = 2.99m~^) 
and for K = —0.04 simulation (short one, Ax = 4.79m~^). 



5. Initial elliptic orbit 



As discussed in Sec. IV B 5[ for e <C 1, the initial orbit of the AD field becomes elliptic in the phase space. The 



resultant extreme energy-to-charge ratio E/mQ ~ 1/e 1 becomes smaller as a result that the Q-balls discard the 
excessive energy. In this process, the negative Q-balls are produced to conserve the total charge. 

We have performed numerical simulations with e = 0.1 (3D4) and e = 0.01 (3D5). The snapshots of the spatial 
charge distribution at r = 1500,2500 and 5000 are shown in Figs. [23] and [24j respectively. The first formation take 
place around r ~ 900, when the filamentary structure is constructed and torn to small pieces, as in the case with 
e = 1. Hence we omit the corresponding plots. The surfaces shown in these figures are the isodensity surfaces with 
q = ±10~^e|<l>inpm, and the red indicates the positive Q-balls and the cyan the negative ones. 

In the panel (a), we can see the first-generation Q-balls. As is same as the 2D cases in Sec. |IVB^ , the all Q-balls 



have positive charges. After that, the energy release from the them is started at r ~ 2000, and as a result, the 
negative Q-balls appear around the first-generation Q-balls, which is shown in the panel (b). This process becomes 
more violent if e is smaller. 

At r = 5000 [the panel (c)], the relaxation process has been almost finished. We can see some isolated negative 
Q-balls, but most of them are still paired with the positive Q-balls. This process to form the second-generation Q-balls 
is quite similar to the cases in the 2D system. Due to the larger spatial dimension, however, it is observed that a 
larger number of the negative Q-balls appeared around one first-generation positive Q-ball. 



Fig. [25] shows the charge distributions of the positive and negative Q-balls. The upper panels are the result with 
6 = 0.1, and the lower ones that with e = 0.01. We take the critical charge Qc = ±1.8 x 10~^e|^inpm. In the 
case of e = 0.1, almost all the Q-balls are positive, and negative Q-balls are subdominant. The peak charge is 
approximately Qp^^k ^ 2 x 10~^|<l>inpm~^ which indicates the scaling as Qp^^k ^ ^- (Recall that the peak charge 
is Qp^ak ^ 0.02|<l>inpm~^ for e = 1. See Fig. 19 ) On the other hand, in the case of e = 0.01, almost same number 



of positive and negative Q-balls are produced after r ~ 2000. The peak charge is approximately 2 x 10 ^ until 
r < 2000 which agrees with the scaling. However, the scaling of the peak charge becomes invalid at r = 5000. We 
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(a) T = 1500 (b) T = 2500 (c) r = 5000 



FIG. 23: The isodensity surface of the charge density with q — ±10~^|$in|^m in the case of e = 0.1. Note that the absolute 
value of this criterion is 60 times larger than qc used in the following analysis for e = 0.1 in order to make the plots clear. The 
red represents the positive Q-balls, and the cyan the negative ones. For these visualizations, we used OpenDX. 




(a) r = 1500 (b) r = 2500 (c) r = 5000 



FIG. 24: The isodensity surface of the charge density with q — ±10~^|$in|^Tn in the case of e = 0.01. Note that the absolute 
value of this criterion is 60 times larger than qc used in the following analysis for e = 0.01 in order to make the plots clear. The 
red represents the positive Q-balls, and the cyan the negative ones. For these visualizations, we used OpenDX. 



have confirmed that it is the Q-balls at the intersections that give dominant contribution to the peak of the charge 
distributions. In particular, large positive and negative Q-balls are found at the intersections in the case of e = 0.01. 

To see the relaxation process where the excited Q-balls change to the states having lower values of E/Q^ we show 
the energy and the charge of the Q-balls in the left panels of Fig. 
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and |27| and the relation between the size 
(diameter) and the charge in the right panels. In the figures, the green symbols ('*') represent the energy and charge 
of the first-generation (excited) Q-balls at r = 1500, implying E/mQ ~ 1/e, and the red symbols ('+') and the blue 
ones ('x') represent the positive and the negative Q-balls at r = 5000, respectively. 

As is similar to the 2D system discussed in Sec. |IVB"5l we can see that the energy-to-charge ratio decreases as time 
goes, while roughly keeping the proportionality between E and Q, and the sizes become almost universal. Looking 
closely, we found that the minimum energy-to-charge ratio is approximately E/mQ 2 at r = 5000, a nd there seems 
to be a Earner which forbids Q-balls to take a smaller value of E/mQ. As is pointed out in Sec. IVB 5, this ratio may 
be determined by the highly non-linear interaction during the secondary formation process. The relaxation process 
in which the value of E/mQ eventually becomes 1 proceeds quite slowly even though there is such a process. Hence 
it can be concluded that, in general, the second-generation Q-balls acquire values of E/mQ ~ 0(1). 

To summarize this section, we have observed that the filamentary structure plays an essential role both in the Q-ball 
formation and in the final charge distribution. Large (positive and negative) Q-balls are created at the intersections 
of the filaments, while small ones in the sides and voids. 
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Q/(l<I>i„lV2) 

FIG. 25: Charge distributions of Q-balls with qc — ±1.8 x 10~^e|$in|^m. The upper panel is the case with e = 0.1 [3D4], and 
the lower one is e = 0.01 [3D5]. The left panels are distributions of negative Q-balls and the right panels are those of positive 
Q-balls. 

V. APPLICATIONS 

In this section we show how the Q-ball charge distribution obtained in our numerical simulations affects the 
cosmological impact of the Q-balls. We consider the decay and evaporation processes in the following. 



A. Fitting formula 



First let us give a fitting formula of the charge distribution. For simplicity we hereafter focus on the case of e = 1. 
The distribution function /Ar(t,Q) is defined such that the number density of Q-balls with a charge between Q and 
Q + at the time t is equal to /Ar(t, Q)SQ. In terms of N{Q), the number of Q-balls with a charge Q in a comoving 
volume, it is written as 



fN{t,Q) = a{t)-^N{Q), 



(9) 



Recall that the scale factor is normalized as a(tin) = 1, where Un = 2/(3m). The numerical result of N{Q)Q is shown 



in Fig. 19 In order to obtain a fitting formula, we adopt the following form. 



a(tYfN(t,Q) = N{Q)=xim^Q 



exp [-xs ( 



(10) 



where Q = (5/(|<l>inp/m^). Assuming that the dominant source for the statistical error comes from N following the 
Poisson distribution, we have fitted the above form to the charge distribution at r = 5000 shown in Fig. [T9j and 
obtained obtained xi 71.2, X2 — 1.29, and xs 1.86 x 10^. The fitting function and the charge distribution are 
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FIG. 26: The charge and the size (diameter d = 2R) of Q-balls with e = 0.1, = ±1-8 x 10"^|$in|^m [3D4]. The green 
represents those of positive Q-balls at r = 1500. The red '+' and the blue 'x' represent positive Q-balls and negative ones at 
T = 5000, respectively. The arrows standing on the horizontal axis represent the grid size at r = 5000 (Ax = 2.99m~^). 




FIG. 27: The charge and the size (diameter d = 2R) of Q-balls with e = 0.01, = 1.8 x 10~^|$in|^m [3D5]. The green 
represents those of positive Q-balls at r = 1500. The red and the blue 'x' represent positive Q-balls and negative ones at 
T = 5000, respectively. The arrows standing on the horizontal axis represent the grid size at r = 5000 (Ax = 2.99m~^). 



shown in Fig. 28 where the vertical axis is a{t)'^ fjsj{t^Q)Q at r = 5000, being equivalent to Fig. 19 Note that the 
above fitting formula is valid only for e = 1 and the charge Q between 0.001 < Q < 0.1. According to the fitting 
formula, the charge distribution a{t)^ fN{t^Q)Q at r = 5000 is peaked at 



Q 



3D,e=l 
peak 



1.9 X 10" 



m 



(11) 



which is greater than the result of Ref. ^ by about 60%. 
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FIG. 28: Charge distribution of Q-balls with parameter set 3D1 and = 1-8 x 10 ^|$in|^^ at r = 5000. The red solid is the 
charge distribution and the green bars represents the Poisson error of the number of Q- balls in each bin, namely, ±V^- The blue 
dotted line is the fitting function given as a^fQ{Q)Q where fN^t^Q) is defined in Eq. (10), and Q denotes Q = (3/(|$in|^/m^). 



Decay 



The decay of the Q-balls proceeds only in the vicinity of the surfaces, and the decay rate is known to be bounded 
above fO^, 



1 
Q 



dQ 



dt 



< 



1927r2 ' 



(12) 



where A is the surface area and u denotes the energy per unit charge of the Q-ball, respectively. According to Ref. , 
the decay rate is indeed almost saturated. In the case of the gravity-mediation type Q-ball, the decay rate is then 
given by 



m 



Q 



(13) 



247t\K\Q' 

87r/{\K\m'^). Note that the decay rate is a decreasing function of Q; 



where we have used uj ^ m and A 
larger Q-balls are more long-lived. 

Now let us estimate the Q-ball contribution to the total energy at the decay. Those Q-balls between Q and Q + 5Q 
decay when H = rQ(Q), emitting the radiation density Spq{Q): 



SpQ _ rnQ fN{TQ^,Q)SQ 



Ptotal 



3r2M2 



(14) 



where ptotai denotes the total energy density at the decay. Using (10), the distribution at the decay is given by 



(15) 



where Hr denotes the Hubble parameter at the reheating. Thus the fractional contribution is proportional to 
Q^/'^N{Q)^ which is multiplied by an additional factor Q^^^ with the charge distribution. Although the peak may 
be shifted due to the additional factor Q^^^, we may neglect it in practice. We may use the charge at the peak to 
estimate the decay rate of the Q-ball. The decay rate then turned out to be about 60% smaller than the existing 
result 
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C. Evaporation 

Through the interaction with the ambient thermal plasma, some amount of charges are evaporated from Q-balls j2TJ 
[22] , For a given thermal history, we can estimate the evaporated charge AQ. For the gravity- mediation type Q-balls, 
the evaporated charge is given by 

if the reheating temperature is higher than the weak scale. 

Before our work, it was customary to assume that almost all the charge is absorbed into the largest Q-balls with a 
charge close to Qp^^k- Then the fraction of the evaporated charge is estimated by AQ/Qp^^k- Since we have obtained 
a fitting formula for the charge distribution, fniQ) (the time variable is dropped for simplicity), we can estimate the 
fraction more precisely as 

/min [Q,^Q]fN{Q)dQ 

jQfNiQ)dQ • ^ ' 



Using the result (10), the fraction of the evaporated charge reads ~ 0.12 and 0.22 for AQ/Qp^^k = 0-1 0-2, 
respectively. Thus the evaporated charge turns out to be O(10)% larger than the previous result.^ Since the fitting 
formula may underestimate the charge distribution for Q < O.lQp^^^j^, the discrepancy may be actually much greater 
for AQ/Q^D ^ < 0.1. If this is the case, it may be an obstacle for realizing a large hierarchy between the baryon and 
lepton asymmetries through the Q-balls [23] , 



VI. CONCLUSION 



In this paper, we have studied physical properties of the gravity- mediation type Q-balls in the AD mechanism on 
ID, 2D and 3D lattice simulations with the cosmic expansion. 

We have found that including cosmic expansion has a crucial effect on the Q-ball formation and its subsequent evo- 
lution, using the ID lattice simulations. In the non-expanding background, the Q-ball interactions through collisions 
never decouple until one large Q-ball is formed. In fact, we confirmed the finite box effect in our simulations without 
cosmic expansion. On the other hand, in the expanding background, the major collisions among large Q-balls do not 
occur frequently. Even in the cases of the 2D and 3D simulations, we expect that, without the cosmic expansion, the 
Q-ball evolution would be quite different. 

In various setups with the cosmic expansion, we have performed simulations in 2D and 3D systems to obtain the 
charge distributions of the Q-balls. For the initial circular orbit with e = 1, the charge at the peak of the distribution 
is Qlf^^ 0, 0.1|$i„|2m-i and Ql^^ 1.9 x 10-^\<^>,^\^m-\ 

We have found that filamentary structure is formed when the initial fluctuation becomes of order unity, as pointed 
out in Ref. [12]. One of our new observations is that this filamentary structure plays an important role in the Q-ball 
formation and affects the final charge distribution. In particular, the Q-balls forming the peak of the distribution 
are mainly produced at the intersections of the filaments, while Q-balls produced from the side and void region have 
generically small charges and do not give main contribution to the charge distribution. Note however that these small 
Q-balls overwhelm the large Q-balls in the number. This argument holds true for a small value of e, in which large 
positive and negative Q-balls are mostly produced from the excited Q-balls at the intersections. (We call the latter 
as the first generation Q-balls and the former as the second generation Q-balls in the text.) 

We have compared the empirical formula derived in Ref. [T2 to our numerical results on the charge (not number) 
distribution on the 2D lattice and found that, if we fit the formula at the peak of the charge distribution, it underes- 
timates the charge distribution at smaller Q. In contrast to the claim made by Ref. [12], therefore, we conclude that 
these small Q-balls do not play any significant role to determine the Q-ball distribution. We also mention here that 
we stopped following the evolution once the Q-balls become no longer resolved compared with our simulation grid 
size, in order to avoid mis-identification of background fluctuations with the Q-balls as well as unphysical effects on 
the Q-ball evolution. This also may be the reason for the difference of ours from the results in the past. 



^ Note that we here focus on the difference arising only from the crude approximation of neglecting the charge distribution. If we take 
account of the difference in Qp^^k' discrepancy becomes larger, but in the opposite direction. 
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Our findings, especially detailed charge distribution of the Q-balls, will be very important for the AD baryogenesis 
and related topics, because the Q-balls are an essential ingredient of the mechanism. For instance, the Q-ball decay 
temperature is modified accordingly, which affects the resultant dark matter density in a scenario that the lightest 
SUSY particles are produced from the Q-ball decay. Also the charge distribution gives a lower bound on the evaporated 
charge, which in turn makes it difficult to realize a large hierarchy between the baryon and lepton asymmetries in the 
scenario using the L-balls. 

Lastly we mention the limitation of our numerical simulations. Strictly speaking, the Q-balls we identified in the 
simulations are in a excited state, in a sense that the E /mQ ratio is slightly larger than 1. We expect that these Q- 
balls will in the end relax to the Q-ball solution, emitting the excessive energy. Also we have neglected the interactions 
of the AD field with gauge fields and other matter fields. We leave these issues as well as further exploration of the 
application of our results for future work. 
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Appendix A: Numerical schemes 



The Q-balls are formed through the dynamical instability of the tiny fluctuations induced at the initial time. 
Until the formation epoch of the Q-balls, it takes a long time, typically a few hundreds of oscillation periods of the 
background field. Hence a higher-order numerical scheme is required so that the results does not suffer from the 
global error when we follow the evolution of the AD field. In this paper, we used a 6th-order symplectic integrator 
developed by Yoshida [16]. 

If we write the time evolution for a time step At as a mapping, ^{t + At^x) = S{At)^{t^x)^ the well-known 
2nd-order leap-frog scheme can be described as 

52nd(At) = e5^*^e^*^e5^*^, (Al) 

where At is the time step, and the operators e^^^^ and e^^^^ map the field and its time-derivative at a time slice t 
to those at the next time slice t + cAt, 

gCAtA . ^(^^^ ^ ^(^^ ^ ^^^^ ^ ^(^^^ ^ cAt^{t, x) , (A2) 

^cAtB . ^^^^ ^ ^ ^^^^ ^ ^^^^ ^ cAt^{t, x). (A3) 



In Eq. (A3), the second derivative of ^ is replaced by i> and ^ using the equation of motion of $ [see Eq. (l2|]. Using 



this symbolic representation, we can describe the 6th-order scheme as [16] 

Seth{At) = S2nd{w3At)S2nd{w2At)S2nd{wiAt)S2nd{woAt)S2nd{wiAt)S2nd{^^ (A4) 



where the coefficients wi,W2 and ws are given in Table |IV| and wq = 1 — 2{wi -\- W2 -\- ws). As implied in Eq. ( |A3[ ), 
the function evaluation, V\^) in the present case, is included only in e^^^^. Therefore this scheme is 7 stages. 



Wi 
W3 



-1.17767998417887 
0.235573213359357 
0.784513610477560 



TABLE IV: A set of coefficients used in 6th-order symplectic scheme I16J. 



As for the spatial derivative, i.e. the Laplacian operator in Eq. ([2]), we used the finite difference scheme called as 
3-, 5- and 7-point formulas for ID, 2D and 3D, respectively. For example, in the 3D cases, we replace the Laplacian 
by 

where ^ij^k ^{t, Xiji/j ^ Zk) ' This derivative is computed at each stage just before the function evaluation. 
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Appendix B: identification of Q-balls 



In order to identify a Q-ball in the simulation box, we developed an identification method based on the Marching 
Cube algorithm with some modifications. The Marching Cube method has been frequently used in imaging analyses 

HZ]. 

We obtain the field data, <l>(t,ic), from the simulations, and can easily compute the charge density, g(t,ic), on 
the grid. In this charge density field, we regard the region in which the charge density exceeds the critical value, 
\q{t,x)\ > as a Q-ball. In accordance with the sign, we call positive Q-ball {q > 0) and negative Q-ball {q < 0). 

Looking closely at a grid cell containing a portion of Q-ball boundaries, the boundary piece lies in the grid cell as 



shown, for example, in Fig. 29 The left cell in the figure is an example of the case where there is only a grid point 
exceeding the critical value, represented by red sphere. The cyan triangle is a piece of the boundary of a Q-ball, 
and the gray region represents the inner region of the Q-ball. The right cell in the figure is an example of a case 
where there are four points exceeding the critical value. The location of the boundary, namely, the spatial coordinate 
satisfying q{t, x) = q^ on the sides of the grid cell, is determined by the linear interpolation between the adjacent grid 
points. 

In accordance with the spatial dimension, there are several patterns of the shape of Q-ball boundary pieces in a 
grid cell. In the ID case, there are only 1 pattern in which the boundary cuts a line segment (ID cell). In the 2D 
cases, we have 4 patterns: a triangle (or, inversely, a pentagon), a trapezoid, a hexagon, and two triangles in a 2D 
cell. Note that the last pattern indicates that there are parts of two Q-balls in the cell. In the 3D cases, we have 14 



patterns including Fig. 29 The all patterns in 3D are shown in Figure 3 in Ref. [171. 

Our algorithm identifying Q-balls is quite simple. First we seek local peaks of the charge density field satisfying 
\q{t^ x)\ > qc^ which are candidates of center of individual Q-balls. Then we seek the boundary outward from the cell 
containing the peak. At the same time, the energy density and charge density of each grid cell are computed. Coming 
to the boundary, we classify the grid cell in the patterns mentioned in the previous paragraph, and compute the 
location of the boundary piece. Finally, integrating the charge density and energy density over the region identified 
as a Q-ball, we obtain the charge, energy and the (D-dimensional) volume of the Q-ball. The size (diameter) of the 
Q-ball is calculated assuming that the shape of the Q-ball is disk (2D) or ball (3D) with the same area (2D) or volume 
(3D). Hence the diameter of the Q-balls is given by 

V for ID, 

d = } ^fWJ^ for 2D, (Bl) 
^^/WJ^ for 3D, 

where V is the D-dimensional volume. 

If by chance there are more than two peaks of the charge density field in the region satisfying \q(t^x)\ > q^^ the 
present scheme counts redundantly multiple Q-balls. Hence, at the final stage of this algorithm, the redundant Q-balls 
are eliminated. 




FIG. 29: Two examples of grid cells containing a piece of boundary of a Q-ball. The left cell is the case where there is only 
a grid point exceeding the criterion in the cell, represented by red sphere, and the right one is a case where there are four 
points exceeding the criterion. The cyan triangle (left) and hexagon (right) are a piece of Q-ball boundary, and the gray region 
represents the inner region of a Q-ball. 
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